#########################
#### Call libraries. ####

library(ssanv)
library(exactci)
library(exact2x2)
library(sensitivitymw)

library(iterators)
library(foreach)
library(doParallel)
library(doMC)
library(plyr)

#######################################
#### Read in the matched datasets. ####

mod_match_ortho_0 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_matched_mod_pats.csv")
mod_match_gensurg_0 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_matched_mod_pats.csv")
both_0 <- rbind.fill(mod_match_ortho_0, mod_match_gensurg_0)

mod_match_ortho_1 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_1_ortho_matched_mod_pats.csv")
mod_match_gensurg_1 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_1_gensurg_matched_mod_pats.csv")
both_1 <- rbind.fill(mod_match_ortho_1, mod_match_gensurg_1)

mod_match_ortho_2 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_ortho_matched_mod_pats_12.21.2018.csv")
mod_match_gensurg_2 <- read.csv("/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_2_gensurg_matched_mod_pats_1.2.2019.csv")
both_2 <- rbind.fill(mod_match_ortho_2, mod_match_gensurg_2)

both_0 <- both_0[,-which(names(both_0) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
both_1 <- both_1[,-which(names(both_1) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
both_2 <- both_2[,-which(names(both_2) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]

#########################################
#### Read in the unmatched datasets. ####

new_ortho_0 <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_new_pats.csv")
exp_ortho_0 <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_ortho_exp_pats.csv")
new_gensurg_0 <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_new_pats.csv")
exp_gensurg_0 <- read.csv(file="/Users/sharpe/Desktop/Young Surgeons/taper_data/lv_0_gensurg_exp_pats.csv")

mod_gensurg_unmatched <- rbind.fill(new_gensurg_0, exp_gensurg_0)
mod_ortho_unmatched <- rbind.fill(new_ortho_0, exp_ortho_0)

mod_gensurg_unmatched <- subset(mod_gensurg_unmatched, era_mod==1)
mod_ortho_unmatched <- subset(mod_ortho_unmatched, era_mod==1)


mod_gensurg_unmatched <- mod_gensurg_unmatched[,-which(names(mod_gensurg_unmatched) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]
mod_ortho_unmatched <- mod_ortho_unmatched[,-which(names(mod_ortho_unmatched) %in% c("icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other"))]

mod_gensurg_unmatched$new_surgeon_flag <- ifelse(mod_gensurg_unmatched$Mod_New_Surgeon_flag_numeric==1,1,0)
mod_ortho_unmatched$new_surgeon_flag <- ifelse(mod_ortho_unmatched$Mod_New_Surgeon_flag_numeric==1,1,0)

###############################################
#### Read in the updated icu_yes variable. ####

trad_ortho_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/trad_ortho_new_cost_and_icu.csv")
mod_ortho_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/mod_ortho_new_cost_and_icu.csv")
trad_gensurg_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/trad_gensurg_new_cost_and_icu.csv")
mod_gensurg_icu <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/mod_gensurg_new_cost_and_icu.csv")

trad_ortho_icu <- trad_ortho_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
mod_ortho_icu <- mod_ortho_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
trad_gensurg_icu <- trad_gensurg_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
mod_gensurg_icu <- mod_gensurg_icu[,c("newID_1","newID_10","icu_yes","total_cost_complete", "total_cost_index_only","icu_days","total_cost_other")]
new_icu_data <- rbind(trad_ortho_icu, mod_ortho_icu, trad_gensurg_icu, mod_gensurg_icu)

both_0 <- merge(x=both_0, y=new_icu_data, by=c("newID_1","newID_10"),all.x=T)
both_1 <- merge(x=both_1, y=new_icu_data, by=c("newID_1","newID_10"),all.x=T)
both_2 <- merge(x=both_2, y=new_icu_data, by=c("newID_1","newID_10"),all.x=T)

mod_gensurg_unmatched <- merge(x=mod_gensurg_unmatched, y=new_icu_data, by=c("newID_1","newID_10"), all.x=T)
mod_ortho_unmatched <- merge(x=mod_ortho_unmatched, y=new_icu_data, by=c("newID_1","newID_10"), all.x=T)

####################################################
##### Read in the readms_or_death_new variable. ####
tg_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_gensurg_new_readms.csv")
mg_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_gensurg_new_readms.csv")
to_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_trad_ortho_new_readms.csv")
mo_readms <- read.csv("/Users/sharpe/Desktop/Young Surgeons/data/transfer_mod_ortho_new_readms.csv")
new_readms_data <- rbind.fill(tg_readms, mg_readms, to_readms, mo_readms)

both_0 <- merge(x=both_0, y=new_readms_data, by=c("newID_1","newID_10"),all.x=T)
both_1 <- merge(x=both_1, y=new_readms_data, by=c("newID_1","newID_10"),all.x=T)
both_2 <- merge(x=both_2, y=new_readms_data, by=c("newID_1","newID_10"),all.x=T)

mod_gensurg_unmatched <- merge(x=mod_gensurg_unmatched, y=new_readms_data, by=c("newID_1","newID_10"), all.x=T)
mod_ortho_unmatched <- merge(x=mod_ortho_unmatched, y=new_readms_data, by=c("newID_1","newID_10"), all.x=T)
mod_unmatched <- rbind.fill(mod_gensurg_unmatched, mod_ortho_unmatched)

#########################
#### Create subsets. ####

gensurg_0 <- subset(both_0, gensurg_adm_mod==1)
gensurg_1 <- subset(both_1, gensurg_adm_mod==1)
gensurg_2 <- subset(both_2, gensurg_adm_mod==1)

ortho_0 <- subset(both_0, ortho_adm_mod==1)
ortho_1 <- subset(both_1, ortho_adm_mod==1)
ortho_2 <- subset(both_2, ortho_adm_mod==1)

######################################################
#### Define the outcomes that we want to look at. ####

bin_outcomes <- c("readms_30","readms_or_death","dead_inhosp","ftr_30_ign_poa","ftr_30_use_poa","ftr_inhosp_ign_poa","ftr_inhosp_use_poa","icu_yes",
                  grep(x=names(both_0),pattern="complic_inhosp_",value=TRUE), "dead_30","dead_60","dead_90","dead_180","prolonged_d","step_up","complic_poa_any",
                  "readms30_post", "readms_or_death_new")
bin_outcomes <- bin_outcomes[-which(bin_outcomes=="complic_inhosp_COMPSYN")]

con_outcomes <- c("total_cost_complete","total_cost_index_only","anesth_time_claim","icu_days","los_long", "basic_los")
# con_outcomes <- c("total_cost_complete", "total_cost_index_only")

####################################################
#### Write function to perform binary outcomes. ####



mcnemar_test <- function(data, treatment, comparison, variable, pair_id){
  
  
  test_data_t1 <- data[which(data[,treatment]==1),c(pair_id,variable)]
  test_data_t0 <- data[which(data[,treatment]==0),c(pair_id,variable)]
  
  colnames(test_data_t1)[2] <- c("var_t1")
  colnames(test_data_t0)[2] <- c("var_t0")
  
  test_data <- merge(x=test_data_t1, y=test_data_t0, by=c(pair_id))
  
  mctable <- table(test_data$var_t0, test_data$var_t1)
  test <- mcnemar.test(x=mctable)
  ci <- mcnemar.exact(x=mctable, conf.level=0.95)
  print(mctable)
  
  #### Get rate of binary variable. ####
  
  overall_rate <- mean(data[,variable], na.rm=TRUE)
  t1_rate <- mean(data[which(data[,treatment]==1),variable],na.rm=TRUE)
  t0_rate <- mean(data[which(data[,treatment]==0),variable],na.rm=TRUE)
  
  mcnemar_test_vector <- c(comparison, variable, test$statistic, test$p.value, test$method, ci$estimate, ci$conf.int, overall_rate, t1_rate, t0_rate)
  names(mcnemar_test_vector) <- c("comparison", "Outcome", "statistic (Chi-square)","McNemar p.value","method","Paired Odds Ratio", "Odds Ratio Lower 95%","Odds Ratio Upper 95%", paste("New and Exp.", variable, "Rate", sep=" "), "New Surgeon Rate", "Exp. Surgeon Rate")
  
  
  return(mcnemar_test_vector)
  
}

#####################################################
#### Write function to perform continuous tests. ####

mestimates <- function(data, treatment, pair_id, variable, comparison){
  
  test_data_t1 <- data[which(data[,treatment]==1),c(pair_id,variable)]
  test_data_t0 <- data[which(data[,treatment]==0),c(pair_id,variable)]
  
  colnames(test_data_t1)[2] <- c("var_t1")
  colnames(test_data_t0)[2] <- c("var_t0")
  
  test_data <- merge(x=test_data_t1, y=test_data_t0, by=c(pair_id))
  test_data <- as.matrix(subset(test_data, select=c("var_t1","var_t0")))
  
  m1 <- senmwCI(y=test_data_t1[,c("var_t1")], gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  m0 <- senmwCI(y=test_data_t0[,c("var_t0")], gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  diff <- senmwCI(y=(test_data_t1[,c("var_t1")] - test_data_t0[,c("var_t0")]), gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  
  if(diff$PointEstimate[1] < 0){
    
    test_data <- test_data*(-1)
    
  }
  
  
  m <- senmw(y=test_data, gamma=1, method="h", inner=0, trim=1, lambda=0.98, tau=0, m=1, m1=1, m2=1)
  
  
  
  m_vec <- c(comparison, variable, (m$pval)*2, m$deviate, m$statistic,m$expectation, m$variance, m1$PointEstimate[1], m1$Confidence.Interval[1],m1$Confidence.Interval[2],m0$PointEstimate[1], m0$Confidence.Interval[1], m0$Confidence.Interval[2], diff$PointEstimate[1], diff$Confidence.Interval[1], diff$Confidence.Interval[2])
  names(m_vec) <- c("Comparison","Outcome","M Statistic P-value", "Deviate", "Statistic", "Expectation", "Variance", "New Surgeons M-estimate", "New Surgeons Lower 95% CI", "New Surgeons Upper 95% CI", "Exp. Surgeons M-estimate", "Exp. Surgeons Lower 95% CI", "Exp. Surgeons Upper 95% CI", "New - Exp. M-estimate",
                    "New - Exp. Lower 95% CI", "New - Exp. Upper 95% CI")
  
  
  
  return(m_vec)
  
  
}

mestimates_overall <- function(data, treatment, variable, comparison){
  
  test_data_t0 <- data[which(data[,treatment]==0),c(variable)]
  test_data_t1 <- data[which(data[,treatment]==1),c(variable)]
  
  
  
  m0 <- senmwCI(y=test_data_t0, gamma=1, method="h", inner=0, trim=1, lambda=0.98, m=1, m1=1, m2=1, alpha=0.05, one.sided=F, tol=NULL, interval=NULL, detail=T)
  m1 <- senmwCI(y=test_data_t1, gamma=1,method="h",inner=0,trim=1,lambda=0.98,m=1,m1=1,m2=1,alpha=0.05,one.sided=F,tol=NULL, interval=NULL, detail=T)
  m_vec <- c(comparison, variable, m0$PointEstimate[1], m0$Confidence.Interval[1],m0$Confidence.Interval[2],  m1$PointEstimate[1], m1$Confidence.Interval[1],m1$Confidence.Interval[2])
  names(m_vec) <- c("Comparison","Outcome","Exp. Surgeon M-Estimate", "Exp Surgeon - Lower 95% CI", "Exp Surgeon - Upper 95% CI","New Surgeon M-Estimate", "New Surgeon - Lower 95% CI", "New Surgeon - Upper 95% CI")
  
  return(m_vec)
}

#######################
#### Do the tests. ####

mcnemar_table <- NULL

for(i in 1:length(bin_outcomes)){
  
  print(bin_outcomes[i])
  mcnemar_both_0 <- mcnemar_both_1 <- mcnemar_both_2 <- mcnemar_gensurg_0 <- mcnemar_gensurg_1 <- mcnemar_gensurg_2 <- mcnemar_ortho_0 <- 
    mcnemar_ortho_1 <- mcnemar_ortho_2 <- NULL
  
  
  mcnemar_both_0 <- mcnemar_test(data=both_0, treatment="new_surgeon_flag", comparison="Level 0 - Gensurg + Ortho", variable=bin_outcomes[i], pair_id="pair_id_pats")
  mcnemar_both_1 <- mcnemar_test(data=both_1, treatment="new_surgeon_flag", comparison="Level 1 - Gensurg + Ortho", variable=bin_outcomes[i], pair_id="pair_id_pats")
  mcnemar_both_2 <- mcnemar_test(data=both_2, treatment="new_surgeon_flag", comparison="Level 2 - Gensurg + Ortho", variable=bin_outcomes[i], pair_id="pair_id_pats")
  
  mcnemar_gensurg_0 <- mcnemar_test(data=gensurg_0, treatment="new_surgeon_flag", comparison="Level 0 - Gensurg", variable=bin_outcomes[i], pair_id="pair_id_pats")
  mcnemar_gensurg_1 <- mcnemar_test(data=gensurg_1, treatment="new_surgeon_flag", comparison="Level 1 - Gensurg", variable=bin_outcomes[i], pair_id="pair_id_pats")
  mcnemar_gensurg_2 <- mcnemar_test(data=gensurg_2, treatment="new_surgeon_flag", comparison="Level 2 - Gensurg", variable=bin_outcomes[i], pair_id="pair_id_pats")
  
  mcnemar_ortho_0 <- mcnemar_test(data=ortho_0, treatment="new_surgeon_flag", comparison="Level 0 - Ortho", variable=bin_outcomes[i], pair_id="pair_id_pats")
  mcnemar_ortho_1 <- mcnemar_test(data=ortho_1, treatment="new_surgeon_flag", comparison="Level 1 - Ortho", variable=bin_outcomes[i], pair_id="pair_id_pats")
  mcnemar_ortho_2 <- mcnemar_test(data=ortho_2, treatment="new_surgeon_flag", comparison="Level 2 - Ortho", variable=bin_outcomes[i], pair_id="pair_id_pats")
  
  mcnemar_table <- rbind(mcnemar_table, mcnemar_both_0, mcnemar_both_1, mcnemar_both_2, mcnemar_gensurg_0, mcnemar_gensurg_1,
                         mcnemar_gensurg_2, mcnemar_ortho_0, mcnemar_ortho_1, mcnemar_ortho_2)
  
}

write.csv(x=mcnemar_table, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/binary_outcomes_1.2.2019_FIXED.csv")


mest_vec_all <- NULL

matched_list <- vector(length=9, mode="list")
matched_list[[1]] <- both_0
matched_list[[2]] <- both_1
matched_list[[3]] <- both_2
matched_list[[4]] <- gensurg_0
matched_list[[5]] <- gensurg_1
matched_list[[6]] <- gensurg_2
matched_list[[7]] <- ortho_0
matched_list[[8]] <- ortho_1
matched_list[[9]] <- ortho_2
names(matched_list) <- c("Level 0 - Gensurg + Ortho", "Level 1 - Gensurg + Ortho", "Level 2 - Gensurg + Ortho",
                         "Level 0 - Gensurg", "Level 1 - Gensurg", "Level 2 - Gensurg",
                         "Level 0 - Ortho", "Level 1 - Ortho", "Level 2 - Ortho")

start <- Sys.time()
registerDoMC(cores=10)

mest_vec_all <- foreach(i=1:length(con_outcomes),.combine='rbind') %:%
  foreach(j=1:length(matched_list),.combine='rbind') %dopar% {
    
    comparison <- names(matched_list)[j]
    
    a_vec <- mestimates(data=matched_list[[j]], treatment="new_surgeon_flag", pair_id="pair_id_pats", variable=con_outcomes[i], comparison=comparison)
    
    
    
  }

end <- Sys.time()
total_time1 <- end-start
total_time1

write.csv(x=mest_vec_all, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/matched_m_estimates_1.2.2019_FIXED.csv")


mest_vec_overall <- NULL

unmatched_list <- vector(length=3,mode="list")
unmatched_list[[1]] <- mod_unmatched
unmatched_list[[2]] <- mod_gensurg_unmatched
unmatched_list[[3]] <- mod_ortho_unmatched




mest_vec_overall <- foreach(i=1:length(con_outcomes),.combine='rbind') %:%
  foreach(j=1:length(unmatched_list),.combine='rbind') %dopar% {
    
    comparison=ifelse(j==1,"Gensurg + Ortho",ifelse(j==2,"Gensurg", "Ortho")) 
    b_vec <- mestimates_overall(data=unmatched_list[[j]], treatment="new_surgeon_flag",variable=con_outcomes[i], comparison=comparison)
    
  }

write.csv(x=mest_vec_overall, file="/Users/sharpe/Desktop/Young Surgeons/taper_table/unmatched_m_estimates_1.2.2019_FIXED.csv")








